﻿extern alias alglibnet2;

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace ClusterAggregation.ClusterAggregators
{
    using Datum;
    using DataSets;
    using Clusterers;
    using GeneralTools;

    /**
     * Volkovich-Yahalom is a special cluster aggregation algorithm.
     * This algorithm uses both the clustering result AND the original dataset to re-evaluate to which cluster each datum belongs.
     * @see ( http://braude-cluster-project.googlecode.com/files/consensus4%20%281%29.pdf )
     */
    public class CVolkovichYahalom:IConsensus
    {
        /** clusterer to use (must be based on a new metric - no new datums outside of given matrix based metric */
        private IClusterer m_clusterer;
     
        /** smoothing factor */
        private double m_lambda;

        /** 
         * simple constructor
         * @param arg (IClusterer) [IN] new clusterer - if null we'll use Spectral Clustering with 0.67 as smoothing factor)
         * @param lambda (double) [IN]  new smoothing factor - default is 1
         */
        public CVolkovichYahalom(IClusterer arg = null, double lambda = 1)
        {
            m_clusterer = arg;
            if (null == m_clusterer)
            {
                //m_clusterer = new CPamImpl();
                m_clusterer = new CSpectralClusterer(null, 0.3, true);
            }

            m_lambda = lambda;
            if (lambda <= 0)
            {
                m_lambda = 1;
            }            
        }

        public CPartition ensemble(CPartition[] partitions, AData[] data, int k)
        {
            int n = data.Length;
            double[,] aggregatedR = new double[n, n];
            double[,] Ri = new double[n, n];

            CMatrixOperators.resetMat(ref aggregatedR, n, n); // initiate R

            foreach (CPartition p in partitions)
            {
                calculateRi(ref Ri, data, n, p);
                CMatrixOperators.addToFirst(ref aggregatedR, Ri, n, n);
            }
            CMatrixOperators.divMat(ref aggregatedR, n, n, partitions.Length);
       
            CMatrixBasedDistance distfunc = new CMatrixBasedDistance(aggregatedR);
           
            CPartition tmp = m_clusterer.cluster(data,distfunc,k);
            tmp.name = "Volkovich-Yahalom";
            tmp.ANMI = tmp.compareTo(partitions);
            return tmp;
        }

        /**
         * helper function to calculate Ri (a d x d matrix that holds the partitioned data's training metric.
         * @param R (double[,]) REF [out] the returning matrix.
         * @param data (AData[]) [IN] the dataset.
         * @param n (int) [IN] the entire dataset's size.
         * @param p (CPartition) [IN] the partition for which we generate Ri.
         */
        public void calculateRi(ref double[,] R, AData[] data, int n, CPartition p)
        {
            int d = data[0].coordinates.Length;
            int C = p.clusters.Count();
            double[,] A = new double[d, d];
            double[,] covMat= new double[d, d];
            double[,] tmp = new double[d, d];
            double[,] Zi = new double[n, C];
            alglibnet2.alglib.matinvreport rep;
            int info;

            CMatrixOperators.resetMat(ref R, n, n);
            CMatrixOperators.resetMat(ref tmp, d, d);
            CMatrixOperators.resetMat(ref A, d, d);
            CMatrixOperators.resetMat(ref Zi, n, C);

            
            for (int j = 0; j < C; j++)
            {

                int clsSize = p.clusters[j].data.Count();
                double tmpSmoother = m_lambda / (m_lambda + clsSize);

                double[,] clsMat = new double[p.clusters[j].data.Count(), d];
                CMatrixOperators.resetMat(ref clsMat, p.clusters[j].data.Count(), d);
                CMatrixOperators.resetMat(ref covMat, d, d);

                // calculate Xi
                AData Xi = new AData(p.clusters[0].data[0].coordinates.Length);
                Xi.zerodize();

                foreach (AData datum in p.clusters[j].data)
                {
                    Xi.accumulate(datum,p.clusters[j].data.Count());
                }

                // for current cluster we have specific Xi
                // calc tmp = Xi Xit
                CMatrixOperators.resetMat(ref tmp, d, d);

                for (int i = 0; i < d; i++)
                {    
                    for (int l = 0; l < d; l++)
                    {
                        tmp[i, l] = Xi.coordinates[i] * Xi.coordinates[l];
                    }
                }

                // calculate covariance matrix
                alglibnet2.alglib.covm(tmp, out covMat);

                // add (tmp = Xi Xit) to A - never forgetting to mult by clsSize!

                for (int i = 0; i < d; i++)
                {
                    for (int l = 0; l < d; l++)
                    {
                        A[i, l] += clsSize * ( covMat[i, l] + tmpSmoother * tmp[i, l]);
                    }
                }


            } // end of calc Ai ^ -1


            //Ai ^ -1
            alglibnet2.alglib.spdmatrixinverse(ref A, out info, out rep);

            tmp = new double[d, n];
            CMatrixOperators.resetMat(ref tmp, d, n);

            // calc Xt
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < d; j++)
                {
                    tmp[j, i] = data[i].coordinates[j];
                }
            }
            //TODO: calc X, then do a mult function for A * B, and use to calc Ri = X * A * Xt

            CMatrixOperators.multMatrices(out covMat, A, tmp, d, d, n);

            tmp = new double[n, d];
            CMatrixOperators.resetMat(ref tmp, n, d);

            // calc X
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < d; j++)
                {
                    tmp[i, j] = data[i].coordinates[j];
                }
            }

            CMatrixOperators.multMatrices(out R, tmp, covMat, n, d, n);

        }
    }
}
